library(tidyverse)
library(lubridate)
library(themebg)
library(plotly)
library(broom)

x <- dirr::get_rds("../data/final")
# data_wt_avg <- read_rds("../data/final/data_wt_avg.rds")

se_line <- list(color = 'rgba(7, 164, 181, 0.05)')
se_fill <- 'rgba(7, 164, 181, 0.2)'
hep_run %>%
    filter(drip.count == 1,
           order %in% c("ACS", "DVT", "Stroke")) %>%
    ggplot(aes(x = run.time, y = med.rate, color = order)) +
    geom_point(shape = 1) +
    # geom_line() +
    geom_smooth() +
    xlim(c(0, 48)) +
    # facet_wrap(~ order, ncol = 1) +
    theme_bg()
hep_run %>%
    filter(drip.count == 1,
           order %in% c("ACS", "DVT", "Stroke")) %>%
    ggplot(aes(x = duration_hypothermia, y = med.rate, color = order)) +
    geom_point(shape = 1) +
    # geom_line() +
    geom_smooth() +
    xlim(c(-24, 72)) +
    theme_bg()
bolus <- hep_bolus_initiation %>%
    distinct(millennium.id) %>%
    mutate(bolus = TRUE)

ptt_run %>%
    left_join(bolus, by = "millennium.id") %>%
    left_join(hep_protocol[c("millennium.id", "order")], by = "millennium.id") %>%
    dmap_at("bolus", ~coalesce(.x, FALSE)) %>%
    filter(order %in% c("ACS", "DVT")) %>%
    ggplot(aes(x = duration_heparin, y = lab.result, color = order)) +
    geom_point(shape = 1) +
    geom_smooth() +
    xlim(c(-12, 96)) +
    theme_bg()
ptt_run %>%
    # mutate(on_heparin = lab.datetime >= heparin.start) %>%
    left_join(bolus, by = "millennium.id") %>%
    dmap_at(c("bolus", "on_heparin"), ~coalesce(.x, FALSE)) %>%
    filter(!is.na(duration_heparin)) %>%
    ggplot(aes(x = duration_hypothermia, y = lab.result, color = bolus)) +
    # geom_point(aes(shape = on_heparin)) +
    geom_vline(xintercept = 0, color = "light grey") +
    geom_point(shape = 1) +
    geom_smooth() +
    scale_shape_manual(values = c(1, 2, 3)) +
    # xlim(c(-12, 96)) +
    scale_x_continuous(breaks = seq(-24, 120, 12), limits = c(-12, 96)) +
    theme_bg()
bolus <- hep_bolus_initiation %>%
    distinct(millennium.id) %>%
    mutate(bolus = TRUE)

df <- ptt_run %>%
    left_join(bolus, by = "millennium.id") %>%
    left_join(hep_protocol[c("millennium.id", "order")], by = "millennium.id") %>%
    dmap_at("bolus", ~coalesce(.x, FALSE)) %>%
    filter(!is.na(duration_heparin),
           duration_hypothermia > -24,
           duration_hypothermia < 72,
           order %in% c("ACS", "DVT", "Stroke")) %>%
    left_join(hep_drip[c("millennium.id", "start.datetime", "stop.datetime")], by = "millennium.id") %>%
    filter(lab.datetime >= start.datetime - hours(4),
           lab.datetime <= stop.datetime + hours(4))

d_acs <- filter(df, order == "ACS")
d_dvt <- filter(df, order == "DVT")
d_cva <- filter(df, order == "Stroke")

m_acs <- loess(lab.result ~ duration_hypothermia, data = df, subset = order == "ACS")
m_dvt <- loess(lab.result ~ duration_hypothermia, data = df, subset = order == "DVT")
m_cva <- loess(lab.result ~ duration_hypothermia, data = df, subset = order == "Stroke")

plot_ly(x = ~duration_hypothermia) %>%
    add_markers(y = ~lab.result, data = df, split = ~order, marker = list(symbol = "circle-open")) %>%
    # add_markers(y = ~lab.result, data = d_acs, marker = list(symbol = "circle-open")) %>%
    # add_markers(y = ~lab.result, data = d_dvt, marker = list(symbol = "circle-open")) %>%
    add_lines(y = ~fitted(m_acs), data = d_acs, name = "ACS", legendgroup = "acs") %>%
    add_lines(y = ~fitted(m_dvt), data = d_dvt, name = "DVT", legendgroup = "dvt") %>%
    add_lines(y = ~fitted(m_cva), data = d_cva, name = "Stroke", legendgroup = "cva") %>%
    add_ribbons(data = augment(m_acs),
                ymin = ~.fitted - 1.96 * .se.fit,
                ymax = ~.fitted + 1.96 * .se.fit,
                line = list(color = 'rgba(7, 164, 181, 0.05)'),
                fillcolor = 'rgba(7, 164, 181, 0.2)',
                showlegend = FALSE,
                legendgroup = "acs") %>%
    add_ribbons(data = augment(m_dvt),
                ymin = ~.fitted - 1.96 * .se.fit,
                ymax = ~.fitted + 1.96 * .se.fit,
                line = list(color = 'rgba(7, 164, 181, 0.05)'),
                fillcolor = 'rgba(7, 164, 181, 0.2)',
                showlegend = FALSE, 
                legendgroup = "dvt") %>%
        add_ribbons(data = augment(m_cva),
                ymin = ~.fitted - 1.96 * .se.fit,
                ymax = ~.fitted + 1.96 * .se.fit,
                line = list(color = 'rgba(7, 164, 181, 0.05)'),
                fillcolor = 'rgba(7, 164, 181, 0.2)',
                showlegend = FALSE, 
                legendgroup = "cva")

PTT change over time from hypothermia initiation by heparin protocol

hep_bolus %>%
    left_join(hep_start, by = "millennium.id") %>%
    left_join(hypothermia_start, by = "millennium.id") %>%
    filter(!is.na(hypothermia_start)) %>%
    mutate(time_heparin = as.numeric(difftime(med.datetime, heparin.start, units = "hours")),
           time_hypothermia = as.numeric(difftime(med.datetime, hypothermia_start, units = "hours"))) %>%
    filter(time_hypothermia >= -24, time_hypothermia <= 72) %>%
    ggplot(aes(x = time_hypothermia, y = med.dose)) +
    geom_point(shape = 1) +
    theme_bg()

Univariate Comparisons

d <- data_wt_avg %>%
    filter(!is.na(temp.time.wt.avg),
           !is.na(ptt.time.wt.avg))

m <- lm(ptt.time.wt.avg ~ temp.time.wt.avg, data = d)

plot_ly(d, x = ~temp.time.wt.avg) %>%
    add_markers(y = ~ptt.time.wt.avg, split = ~group, marker = list(symbol = "circle-open")) %>%
    add_lines(y = ~fitted(m), showlegend = FALSE) %>%
    add_ribbons(data = augment(m),
                ymin = ~.fitted - 1.96 * .se.fit,
                ymax = ~.fitted + 1.96 * .se.fit,
                line = list(color = 'rgba(7, 164, 181, 0.05)'),
                fillcolor = 'rgba(7, 164, 181, 0.2)',
                showlegend = FALSE)

Time-Weighted Average of Temperature vs. PTT

mod <- lm(ptt.time.wt.avg ~ temp.time.wt.avg, data_wt_avg) 
    
mod %>%
    glance() %>%
    knitr::kable(caption = "Linear Regression Model for Time-Weighted Average PTT by Temperature")
Linear Regression Model for Time-Weighted Average PTT by Temperature
r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC deviance df.residual
0.0026356 -0.0029053 25.27819 0.4756583 0.4912837 2 -845.0907 1696.181 1705.793 115017.6 180
mod %>%
    tidy() %>%
    knitr::kable(caption = "Coefficients for PTT by Temperature")
Coefficients for PTT by Temperature
term estimate std.error statistic p.value
(Intercept) 155.3361788 134.696583 1.1532303 0.2503450
temp.time.wt.avg -0.9894958 1.434718 -0.6896798 0.4912837
d <- temp_ptt_full %>%
    filter(!is.na(vital.result),
           !is.na(lab.result),
           vital.result > 85)

d_study <- filter(d, group == "Study")
d_contr <- filter(d, group == "Control")

# m <- lm(lab.result ~ vital.result, data = d)
m_study <- lm(lab.result ~ vital.result, data = d, subset = group == "Study")
m_contr <- lm(lab.result ~ vital.result, data = d, subset = group == "Control")

plot_ly(x = ~vital.result) %>%
    add_markers(data = d, 
                y = ~lab.result, 
                split = ~group, 
                marker = list(symbol = "circle-open")) %>%
    add_lines(y = ~fitted(m_study), 
              data = d_study, 
              legendgroup = "study", 
              name = "Study") %>%
    add_lines(y = ~fitted(m_contr), 
              data = d_contr, 
              legendgroup = "control", 
              name = "Control") %>%
    add_ribbons(data = augment(m_study),
                ymin = ~.fitted - 1.96 * .se.fit,
                ymax = ~.fitted + 1.96 * .se.fit,
                line = se_line,
                fillcolor = se_fill,
                legendgroup = "study",
                showlegend = FALSE) %>%
    add_ribbons(data = augment(m_contr),
                ymin = ~.fitted - 1.96 * .se.fit,
                ymax = ~.fitted + 1.96 * .se.fit,
                line = se_line,
                fillcolor = se_fill,
                legendgroup = "control",
                showlegend = FALSE)

Temperature vs. PTT

mod <- lm(lab.result ~ vital.result, temp_ptt_full)

mod %>%
    glance() %>%
    knitr::kable(caption = "Linear Regression Model for PTT by Temperature")
Linear Regression Model for PTT by Temperature
r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC deviance df.residual
0.0408339 0.0406413 45.70651 212.0524 0 2 -26115.8 52237.59 52257.13 10405733 4981
mod %>%
    tidy() %>%
    knitr::kable(caption = "Coefficients for PTT by Temperature")
Coefficients for PTT by Temperature
term estimate std.error statistic p.value
(Intercept) 352.027639 19.1660628 18.36724 0
vital.result -2.945295 0.2022587 -14.56202 0
d <- temp_ptt_full %>%
    filter(!is.na(vital.result),
           !is.na(med.rate),
           vital.result > 85)

d_study <- filter(d, group == "Study")
d_contr <- filter(d, group == "Control")

m_study <- lm(med.rate ~ vital.result, data = d, subset = group == "Study")
m_contr <- lm(med.rate ~ vital.result, data = d, subset = group == "Control")

plot_ly(x = ~vital.result) %>%
    add_markers(data = d, 
                y = ~med.rate, 
                split = ~group, 
                marker = list(symbol = "circle-open")) %>%
    add_lines(y = ~fitted(m_study), 
              data = d_study, 
              legendgroup = "study", 
              name = "Study") %>%
    add_lines(y = ~fitted(m_contr), 
              data = d_contr, 
              legendgroup = "control", 
              name = "Control") %>%
    add_ribbons(data = augment(m_study),
                ymin = ~.fitted - 1.96 * .se.fit,
                ymax = ~.fitted + 1.96 * .se.fit,
                line = se_line,
                fillcolor = se_fill,
                legendgroup = "study",
                showlegend = FALSE) %>%
    add_ribbons(data = augment(m_contr),
                ymin = ~.fitted - 1.96 * .se.fit,
                ymax = ~.fitted + 1.96 * .se.fit,
                line = se_line,
                fillcolor = se_fill,
                legendgroup = "control",
                showlegend = FALSE)

Temperature vs. Heparin

temp_ptt_full %>%
    ungroup() %>%
    select(-study_id) %>%
    select_if(is.numeric) %>%
    cor(use = "pairwise.complete.obs") %>%
    tidy() %>%
    knitr::kable(caption = "Correlation Analysis")
Correlation Analysis
.rownames vital.result lab.result med.rate
vital.result 1.0000000 -0.2020739 0.2304837
lab.result -0.2020739 1.0000000 -0.0829357
med.rate 0.2304837 -0.0829357 1.0000000
d <- data_wt_avg %>%
    filter(!is.na(hep.time.wt.avg),
           !is.na(ptt.time.wt.avg))

m <- lm(ptt.time.wt.avg ~ hep.time.wt.avg, data = d)

plot_ly(d, x = ~hep.time.wt.avg) %>%
    add_markers(y = ~ptt.time.wt.avg, split = ~group, marker = list(symbol = "circle-open")) %>%
    add_lines(y = ~fitted(m), showlegend = FALSE) %>%
    add_ribbons(data = augment(m),
                ymin = ~.fitted - 1.96 * .se.fit,
                ymax = ~.fitted + 1.96 * .se.fit,
                line = list(color = 'rgba(7, 164, 181, 0.05)'),
                fillcolor = 'rgba(7, 164, 181, 0.2)',
                showlegend = FALSE)

Time-Weighted Average of Heparin Dose vs. PTT

mod <- lm(ptt.time.wt.avg ~ hep.time.wt.avg, data_wt_avg) 
    
mod %>%
    glance() %>%
    knitr::kable(caption = "Linear Regression Model for Time-Weighted Average PTT by Heparin Dose")
Linear Regression Model for Time-Weighted Average PTT by Heparin Dose
r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC deviance df.residual
0.1044035 0.0970625 24.96059 14.22206 0.0002519 2 -573.8852 1153.77 1162.231 76009.82 122
mod %>%
    tidy() %>%
    knitr::kable(caption = "Coefficients for PTT by Heparin Dose")
Coefficients for PTT by Heparin Dose
term estimate std.error statistic p.value
(Intercept) 91.124251 5.932989 15.358912 0.0000000
hep.time.wt.avg -2.062553 0.546920 -3.771215 0.0002519
d <- temp_ptt %>%
    filter(!is.na(ptt),
           !is.na(hep),
           temp > 85)

# m <- loess(ptt ~ hep, data = d)
m <- lm(ptt ~ hep, data = d)

plot_ly(d, x = ~hep) %>%
    add_markers(y = ~ptt, split = ~group, marker = list(symbol = "circle-open")) %>%
    add_lines(y = ~fitted(m), showlegend = FALSE) %>%
    add_ribbons(data = augment(m),
                ymin = ~.fitted - 1.96 * .se.fit,
                ymax = ~.fitted + 1.96 * .se.fit,
                line = list(color = 'rgba(7, 164, 181, 0.05)'),
                fillcolor = 'rgba(7, 164, 181, 0.2)',
                showlegend = FALSE)

Heparin Dose vs. PTT

mod <- lm(ptt ~ hep, temp_ptt)

mod %>%
    glance() %>%
    knitr::kable(caption = "Linear Regression Model for PTT by Heparin Dose")
Linear Regression Model for PTT by Heparin Dose
r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC deviance df.residual
0.0531423 0.0504675 49.42536 19.86821 1.12e-05 2 -1892.704 3791.409 3803.034 864774.7 354
mod %>%
    tidy() %>%
    knitr::kable(caption = "Coefficients for PTT by Heparin Dose")
Coefficients for PTT by Heparin Dose
term estimate std.error statistic p.value
(Intercept) 60.119281 6.2267289 9.655034 0.00e+00
hep 2.699693 0.6056687 4.457377 1.12e-05
d <- data_wt_avg %>%
    filter(!is.na(hep.time.wt.avg),
           !is.na(temp.time.wt.avg))

m <- lm(hep.time.wt.avg ~ temp.time.wt.avg, data = d)

plot_ly(d, x = ~temp.time.wt.avg) %>%
    add_markers(y = ~hep.time.wt.avg, split = ~group, marker = list(symbol = "circle-open")) %>%
    add_lines(y = ~fitted(m), showlegend = FALSE) %>%
    add_ribbons(data = augment(m),
                ymin = ~.fitted - 1.96 * .se.fit,
                ymax = ~.fitted + 1.96 * .se.fit,
                line = list(color = 'rgba(7, 164, 181, 0.05)'),
                fillcolor = 'rgba(7, 164, 181, 0.2)',
                showlegend = FALSE)

Time-Weighted Average of Temperature vs. Heparin Dose

mod <- lm(hep.time.wt.avg ~ temp.time.wt.avg, data_wt_avg) 
    
mod %>%
    glance() %>%
    knitr::kable(caption = "Linear Regression Model for Time-Weighted Average Heparin Dose by Temperature")
Linear Regression Model for Time-Weighted Average Heparin Dose by Temperature
r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC deviance df.residual
0.0661521 0.0584344 4.000886 8.571428 0.0040808 2 -344.0627 694.1253 702.5619 1936.857 121
mod %>%
    tidy() %>%
    knitr::kable(caption = "Coefficients for Heparin Dose by Temperature")
Coefficients for Heparin Dose by Temperature
term estimate std.error statistic p.value
(Intercept) -62.0534314 24.6203929 -2.520408 0.0130236
temp.time.wt.avg 0.7674734 0.2621421 2.927700 0.0040808
d <- temp_ptt %>%
    filter(!is.na(temp),
           !is.na(hep),
           temp > 85)

m <- lm(hep ~ temp, data = d)

plot_ly(d, x = ~temp) %>%
    add_markers(y = ~hep, split = ~group, marker = list(symbol = "circle-open")) %>%
    add_lines(y = ~fitted(m), showlegend = FALSE) %>%
    add_ribbons(data = augment(m),
                ymin = ~.fitted - 1.96 * .se.fit,
                ymax = ~.fitted + 1.96 * .se.fit,
                line = list(color = 'rgba(7, 164, 181, 0.05)'),
                fillcolor = 'rgba(7, 164, 181, 0.2)',
                showlegend = FALSE)

Temperature vs. Heparin Dose

mod <- lm(hep ~ temp, temp_ptt)

mod %>%
    glance() %>%
    knitr::kable(caption = "Linear Regression Model for Heparin Dose by Temperature")
Linear Regression Model for Heparin Dose by Temperature
r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC deviance df.residual
0.0144323 0.0116482 4.305825 5.183851 0.023394 2 -1023.888 2053.776 2065.401 6563.205 354
mod %>%
    tidy() %>%
    knitr::kable(caption = "Coefficients for Heparin Dose by Temperature")
Coefficients for Heparin Dose by Temperature
term estimate std.error statistic p.value
(Intercept) 26.3372936 7.4747202 3.523516 0.0004817
temp -0.1814172 0.0796805 -2.276807 0.0233940

Multiple Variable Models

mod <- lm(ptt.time.wt.avg ~ temp.time.wt.avg + hep.time.wt.avg, data_wt_avg) 
    
mod %>%
    glance() %>%
    knitr::kable(caption = "Linear Regression Model for Time-Weighted Average PTT by Temperature and Heparin Dose")
Linear Regression Model for Time-Weighted Average PTT by Temperature and Heparin Dose
r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC deviance df.residual
0.1204563 0.1057972 24.31333 8.217189 0.0004523 3 -565.5069 1139.014 1150.262 70936.54 120
mod %>%
    tidy() %>%
    knitr::kable(caption = "Coefficients for PTT by Temperature and Heparin Dose")
Coefficients for PTT by Temperature and Heparin Dose
term estimate std.error statistic p.value
(Intercept) 91.0870879 153.4949841 0.5934206 0.5540166
temp.time.wt.avg 0.0051671 1.6484922 0.0031344 0.9975043
hep.time.wt.avg -2.1647088 0.5524533 -3.9183564 0.0001487
d <- temp_ptt_full %>%
    left_join(hep_bolus_initiation[c("millennium.id", "med.dose")], by = "millennium.id") %>%
    mutate(bolus = !is.na(med.dose))

mod <- lm(lab.result ~ vital.result + med.rate + order + bolus, d)

mod %>%
    glance() %>%
    knitr::kable(caption = "Linear Regression Model for PTT by Temperature and Heparin Dose")
Linear Regression Model for PTT by Temperature and Heparin Dose
r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC deviance df.residual
0.1624535 0.1610094 45.05748 112.4989 0 7 -18222.62 36461.24 36510.49 7065015 3480
mod %>%
    tidy() %>%
    knitr::kable(caption = "Coefficients for PTT by Temperature and Heparin Dose")
Coefficients for PTT by Temperature and Heparin Dose
term estimate std.error statistic p.value
(Intercept) 519.9876272 23.4672951 22.1579703 0.0000000
vital.result -4.6888899 0.2493342 -18.8056412 0.0000000
med.rate -0.1214811 0.1205817 -1.0074591 0.3137842
orderDVT 28.9335230 1.9749763 14.6500606 0.0000000
orderFixed 6.1052397 8.0759564 0.7559773 0.4497140
orderStroke -8.9218505 3.2037181 -2.7848426 0.0053843
bolusTRUE 11.7471425 1.6950694 6.9301837 0.0000000

Chelsea’s Data

plot_ly(ck_heparin, x = ~as.numeric(duration), y = ~Value, split = ~Control)
plot_ly(ck_ptt, x = ~duration_ptt, y = ~Value, split = ~Control) %>%
    add_markers(marker = list(symbol = "circle-open"))
d <- ck_data %>%
    filter(!is.na(Temperature),
           !is.na(heparin))
           # heparin > 0,
           # heparin < 4000)

m <- lm(heparin ~ Temperature, data = d)

plot_ly(d, x = ~Temperature) %>%
    add_markers(y = ~heparin, split = ~Control, marker = list(symbol = "circle-open")) %>%
    add_lines(y = ~fitted(m), showlegend = FALSE) %>%
    add_ribbons(data = augment(m),
                ymin = ~.fitted - 1.96 * .se.fit,
                ymax = ~.fitted + 1.96 * .se.fit,
                line = list(color = 'rgba(7, 164, 181, 0.05)'),
                fillcolor = 'rgba(7, 164, 181, 0.2)',
                showlegend = FALSE)

Temperature vs. Heparin Dose

d <- ck_data %>%
    filter(!is.na(Temperature),
           !is.na(PTT))

d_study <- filter(d, Control == "Study")
d_ctrl <- filter(d, Control == "Control")

m_study <- lm(PTT ~ Temperature, data = d_study)
m_ctrl <- lm(PTT ~ Temperature, data = d_ctrl)

plot_ly(x = ~Temperature) %>%
    add_markers(data = d, y = ~PTT, split = ~Control, marker = list(symbol = "circle-open")) %>%
    add_lines(data = d_study, y = ~fitted(m_study), name = "Study", legendgroup = "study") %>%
    add_lines(data = d_ctrl, y = ~fitted(m_ctrl), name = "Control", legendgroup = "control") %>%
    add_ribbons(data = augment(m_study),
                ymin = ~.fitted - 1.96 * .se.fit,
                ymax = ~.fitted + 1.96 * .se.fit,
                line = list(color = 'rgba(7, 164, 181, 0.05)'),
                fillcolor = 'rgba(7, 164, 181, 0.2)',
                legendgroup = "study",
                showlegend = FALSE) %>%
    add_ribbons(data = augment(m_ctrl),
                ymin = ~.fitted - 1.96 * .se.fit,
                ymax = ~.fitted + 1.96 * .se.fit,
                line = list(color = 'rgba(7, 164, 181, 0.05)'),
                fillcolor = 'rgba(7, 164, 181, 0.2)',
                legendgroup = "control",
                showlegend = FALSE) 

Temperature vs. PTT